Electronic Coherence Dephasing in Excitonic Molecular Complexes: Role of Markov 

and Secular Approximations 
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We compare four different types of equations of motion for reduced density matrix of a system 
of molecular excitons interacting with thermodynamic bath. All four equations are of second order 
in the linear system-bath interaction Hamiltonian, with different approximations applied in their 
derivation. In particular we compare time-nonlocal equations obtained from so-called Nakajima- 
Zwanzig identity and the time-local equations resulting from the partial ordering prescription of the 
cummulant expansion. In each of these equations we alternatively apply secular approximation to 
decouple population and coherence dynamics from each other. We focus on the dynamics of intra- 
band electronic coherences of the excitonic system which can be traced by coherent two-dimensional 
spectroscopy. We discuss the applicability of the four relaxation theories to simulations of popula- 
tion and coherence dynamics, and identify features of the two-dimensional coherent spectrum that 
allow us to distinguish time-nonlocal effects. 
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INTRODUCTION 



Modeling molecular properties related to their non- 
equilibrium dynamics requires various theoretical ap- 
proaches depending on the particular microscopic pro- 
cesses related to the observed molecular features. Since 
the dawn of quantum mechanics, properties of molecules 
and solids have been studied theoretically in ever greater 
detail. This has led in recent years to a state in which 
dynamics of complex systems with multitude of degrees 
of freedom (DOF) is accessible to quantitative theoret- 
ical study Many properties of molecular systems 
are directly related to the equilibrium or time depen- 
dent conformations of nuclear DOF for which electronic 
states play the role of a background contributing to the 
nuclear potential energy surfaces. Problems like these 
are the realm of molecular dynamics (MD) in its clas- 
sical, quantum or mixed versions and quantum chem- 
istry (QC), where impressive qualitative and quantita- 
tive results have been achieved in recent years. For cer- 
tain types of dynamical problems, however, less expen- 
sive model approaches are the preferred choice due to the 
scale of studied system or due to the physical nature of 
studied processes. A good example of such a problem is 
ultrafast photo-induced excited state dynamics of small 
molecular systems and their aggregates [2]. Here, most 
of the relevant experimental information is only avail- 
able through ultrafast non-linear spectroscopy, and thus 
the theory has to span the whole distance between the 
microscopic dynamics of the molecular system, and the 
macroscopic description of experimental signals Q . Typ- 
ical field in which such an approach has yielded deep un- 
derstanding of the relevant physico-chemical processes is 
the study of primary processes in photosynthesis. The re- 
lated quantum mechanical problem is usually formulated 



in terms of a model describing the relevant DOF of the 
system (electronic states of photosynthetic molecules), 
and a thermodynamics bath (the protein environment). 
Parameters for such models can be supplied by experi- 
ment Q , QC studies [H, @ , MD modeling Q , or are a 
result of suitable simplified models [3|. 

Recent advances in non-linear spectroscopy have 
opened a wide new experimental window into the details 
of ultrafast photo-induced dynamics of molecular sys- 
tems. Experimental realization of two-dimensional (2D) 
coherent spectroscopy in the visible and near IR regions 
has enabled to overcome some of the frequency- vs. 
time-resolution competition problem otherwise faced by 
ultrafast spectroscopy, and yielded thus unprecedented 
experimental details of the time evolution of molecular 
excitations. Most importantly, it was predicted that the 
presence of certain oscillatory features in 2D spectra is 
a manifestation of coherences between molecular excited 
states (l2l Il3j . It was also concluded that these oscil- 
lations should be present in the 2D spectrum of photo- 
synthetic Fenna-Matthews-Olson (FMO) chromophore- 
protein complex [12] . Experimental results not only con- 
firmed this prediction [14| . but yielded also surprising 
results such as unexpectedly long life time of these co- 
herences, as compared to the predictions of standard de- 
phasing rate theory. Furthermore, while possible coher- 
ence transfer was ignored by the relaxation theory used in 
Ref. fl2l |. the experiment provided some evidence for its 
role in excitation energy transfer. It was speculated that 
photosynthetic systems might use the coherent mode of 
energy transfer to more efficiently channel excitation en- 
ergy by scanning their energetic landscape in a process 
similar to quantum computing [l4j . More experiments 
have recently reported coherent dynamics in photosyn- 
thetic systems [15] and conjugated polymers [16] . and 
the field of energy transfer in photosynthesis has seen an 
increased interest from theoretical researchers from pre- 
viously unrelated fields (l7l - [20| . 
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Theoretical basis for the description of the decoher- 
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ence phenomena in excitation energy transfer has been 
developed long ago in the framework of the reduced den- 
sity matrix (RDM) pll . l22l | . Equations of motion result- 
ing from this scheme are characterized by the presence 
of time retarded terms responsible for energy relaxation 
and decoherence processes. Equations of this type will be 
denoted as time non-local in this paper. Later, an alter- 
native approach to the derivation of the RDM equations 
of motion has emerged which yields time local equation 
of motion [23l . |24| . Both theories express the relaxation 
term in form of an infinite series in terms of the system- 
bath interaction Hamiltonian, but differ in time ordering 
prescriptions for the cumulant expansion of the evolu- 
tion operator. The time local equations correspond to 
so-called partial time ordering prescription of the cum- 
mulant expansion, while the time non-local equations re- 
sult from so-called chronological time ordering (25l |2(| . 
Although the two schemes yield formally different equa- 
tions of motion for the RDM, they are in fact equivalent 
as long as the complete summation of the corresponding 
infinite series is performed. When the infinite series are 
truncated at a finite order, the two theories yield equa- 
tions that predict different RDM dynamics. This is a 
result of different statistical assumptions about the bath 
that are implicitly made in the two cases (25l . l2a| . In 
all orders of expansion, so-called Markov approximation 
can be used to transform the time non-local equation of 
motion into a certain time-local form. This has to be re- 
garded as an additional approximation which simplifies 
the numerical treatment of the time non-local equations. 
Interestingly, in the second order the time-local equations 
and the time non-local equations with Markov approxi- 
mation have exactly the same form. 

Until recently, most experiments were not sensitive to 
coherence between electronic levels. This allowed a host 
of further approximations to simplify equations of mo- 
tion. Most notably, the secular approximation, which 
amounts to decoupling RDM elements oscillating on dif- 
ferent frequencies from each other, has limited the energy 
transfer phenomena to separate dynamics of population 
transfer and coherence dephasing [27| . Even on very 
short times scale, experiments aimed at studying pop- 
ulation dynamics (pump probe) are not sensitive enough 
to coherences between electronic levels to require non- 
secular theory, although it was suggested that measured 
relaxation time can be distorted by non-secular effects 
(28| . Consequently, most of the theory developed for eval- 
uation of experiments has been aimed at improving calcu- 
lation of the relaxation rates (29l-[3l|. With experiments 
now uncovering new details about the role of electronic 
coherence, theoretical methods beyond rate equations for 
probabilities which are both accurate and numerically 
tractable are required. Although schemes for construct- 
ing equations of motion for the RDM beyond second or- 
der, based on co-propagation of the RDM with auxiliary 
operators, seem feasible and promising [32l |33| . second 
order theories might still be the only option for treatment 
of extended molecular systems. It was suggested previ- 



ously that second order perturbation theory with respect 
to system-bath coupling provides a suitable framework 
for development of such methods [34j ■ This notion is also 
supported by the fact that in the special case of so-called 
spin-boson model, second order time-local equation of 
motion already represents an exact equation of motion 
for the RDM (3^. 

In this paper we study the following four different sec- 
ond order theories: (a) full time non-local (full TNL) 
equation of motion resulting from the Nakajima-Zwanzig 
identity or equivalently from the chronological ordering 
prescription in the cummulant expansion, (b) the full 
time local (full TL) equation of motion resulting from 
the partial ordering prescription in the cummulant ex- 
pansion, or equivalently from Markov approximation ap- 
plied to TNL equation, (c) time non-local equation with 
secular approximation (secular TNL), and (d) time local 
equation with secular approximation (secular TL). We 
discuss the applicability of these equations to the descrip- 
tion of the energy relaxation and decoherence dynamics 
in small systems of molecular excitons with the empha- 
sis of on recent 2D spectroscopic experiments and the 
dynamics of coherence between electronic excited states. 
Note that, in this paper, full refers to the equations where 
no secular approximation has been applied. These equa- 
tions are still of second order of perturbation theory. 

The paper is organized as follows. The next section 
introduces Hamiltonian description of an aggregate of 
small molecules embedded in a protein or solid state 
environment. In Section Mil we describe the details of 
second order theory of system-bath interaction and we 
derive four different equations of motion for the RDM de- 
scribing electronic states of a molecular aggregate. Two- 
dimensional coherent spectroscopy, and non-linear spec- 
troscopy in general are introduced in Section IIV1 In Sec- 
tion |V] we present and discuss numerical results compar- 
ing different theories of relaxation on calculations of co- 
herence life time and 2D spectra. 



II. MODEL HAMILTONIAN 

The investigated molecular system is an aggregate 
composed of N monomers embedded in protein environ- 
ment. Let us first consider a monomeric molecule (a 
chromophore) embedded in the environment, but insu- 
lated from interaction with its neighboring monomers. 
The monomer Hamiltonian has a form 

H m = (4 m) + T(P m ) + Vj m \Q m )) \g m )(g m \ 



+ (4™) +T(P m ) + K M (Q™)) |e m )(e m | , (I) 

where \g m ), |e m ) denote electronic ground and excited 
states, and £g m \ represent electronic energies of 
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these states. The kinetic term T(P m ) and the poten- 
tial terms Vg m \Q m ), Ve m \Q m ) represent the intra- 
molecular DOF and the protein environment (bath or 
reservoir) interacting with these states. By Q m (P m ) we 
denote the (possibly macroscopic) set of coordinates (im- 
pulses) describing both the intramolecular nuclear DOF 
of the m-th monomer as well as the DOF of its surround- 
ings. The total Hamiltonian of the monomer can be split 
into the system, reservoir and the system-reservoir cou- 
pling terms 

H^ = e^\g m )(g m \ 
+ (4 m) + (V e {m) (Q m ) - V^\Q m )) eq )\e m )(e m \ , (2) 
H^^[T{P m ) + V^\Q m )] 

x (\9m){g m \ + |e m )(e m |) 
= [T{P m ) + V^ m \Q m )]®l, (3) 

H^ K ^{v^\Q m )-V^\Q m ) 

-(V e im) (Q m )-V} m HQ m )) eq )\e m )(e m \ 

= A$ m (Q m )|e m )(e m . (4) 

Here, (A(Q)) eq represents averaging of an arbitrary 
Q— dependent operator over equilibrium state of the 
bath. By this choice of the splitting we have assured 
that A$( m )(Q m ) = for the syste m in equilibrium. To 
simplify the notation, we redefine electronic energy of 
the excited state to include the equilibrium average of 
the potential energy difference between the electronic ex- 
cited and ground states, ei" 1 "* = Ee + (Ve m \Qm) — 
Vg m \Q m )) eq and we drop the tilde over ei" 1 ' further in 
this paper. 

An aggregate built out of these monomers can be rep- 
resented on a Hilbert space composed of collective aggre- 
gate states. We define the aggregate ground state 

JV 

\g) = II ®l3m>> ( 5 ) 

m— 1 

states with a single excitation 

n-1 N 

\ u n) = \\ ®\9m) ® |e n ) Yl <8>|sw), (6) 

m— 1 rn'— ra+1 

and multi-excited states in an analogical manner. We 
drop the sign <g> in further consideration for the sake of 



brevity. The Hamiltonian of the aggregate is constructed 
using the energies of collective states 

Trnon — int _ I \ / 

+ ^2(Ae n + n)\u n )(u„\+h. e.t., (7) 

n 

where e g = J2 n £ s™ ) > ^ = N_1 J2 n £ ^ > and = £ ^ ~ 
^ + Em/n4"''- The abreviation h. e. t. denotes higher 
excitonic terms. Due to the fact that the monomers are 
positioned in a tight aggregate, we have to account for 
the interaction energy between their excited states. The 
interaction energy between states \u m ) and \u n ) will be 
denoted J m n, and the corresponding contribution to the 
total Hamiltonian reads 

H s nt = ( J n m Wn)(u m \ + c.c.) + h. e. t.. (8) 

Due to the off-diagonal terms J mn the collective states 
defined in Eq. ((SJ are not eigenstates of the total Hamil- 
tonian H s = R n s on - mt +H s nt . Although the basis of the 
states \g) 1 \u n ) and multiple-excitation states of the ag- 
gregate provides efficient means for defining the Hamil- 
tonian, it is more practical to switch into the basis of 
eigenstates of the Hamiltonian Hs- One of the reasons 
is that while matter interacts with light, the differences 
between eigenenergies of Hs define the resonant transi- 
tion frequencies. The system-bath coupling part of the 
aggregate Hamiltonian reads 

H S -b =^A$ n |u n >(u n |+h. c.t., ( 9 ) 

n 

and thus no terms in the total Hamiltonian couple the 
ground state with the first excited state or higher ex- 
cited state bands. In fact, the total Hamiltonian splits 
into blocks separated approximately by the energy ft£l 
(see Fig. [IJ. This reflects the neglecting of all adiabatic 
couplings, which are supposed to be so weak that they 
do not lead to transitions on the time scale of interest 
(femto and picoseconds). This property is well justified 
e.g. for chlorophyll systems. 

For the subsequent use in this paper, we denote the 
eigenstates of the total electronic Hamiltonian Hs by 
\u a ), a = 1,...,N, for single exciton states formed as 
linear combinations of single excitation states \u n ) and 
\U a ), a = N + 1, . . . , N + N{N - I )/2, for two-exciton 
states formed from the linear combination of pairs of sin- 
gle excitation states. 

III. SECOND ORDER RELAXATION 
THEORIES 

In this section we consider interaction of the electronic 
system described by the Hamiltonian Hs with a macro- 
scopic bath composed of the DOF of the molecular sur- 
roundings. 
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Figure 1: Illustration of the level structure of an excitonic 
system. The excited states \e n ) of N monomers with tran- 
sition frequency Q. (left part of the figure) split due to the 
resonance interaction into N one exciton states \u n ) (right). 
Absorption spectra of ensembles of non-interacting (left) and 
interaction (right) monomers. The system also exhibits higher 
excitons states (two-exciton states \U n ) are depicted here), 
with SI being the mean transition frequency from the one- to 
two-exciton bands. A pictorial 2D spectrum with peaks re- 
sulting from transitions between the ground- and one-exciton 
states (red arrows) and one- and two-exciton states (blue ar- 
rows) is presented in the upper left corner of the figure. The 
transitions between the ground- and one-excitons states lead 
to positive contributions to the 2D spectrum (absorption and 
ground state bleach), while the transitions between the one- 
and two-exciton states result in a negative contribution (ex- 
cited state absorption). 



Projection operator technique and 
Nakajima-Zwanzig identity 



the system, the bath and their interaction will be used 
later in this paper. A formal solution of Eq. (fTT)|) can be 
written using evolution superoperators such that 



W{t)=U{t-t )W(t ) 



(11) 



We are often not interested in the complete information 
carried by W(t), but only in the information related to 
the electronic states of the aggregate. The reduced den- 
sity matrix (RDM) defined as 



p(t)=tr B {W(t)}, 



(12) 



where trs represents trace over all bath DOF, is the 
quantity carrying just the relevant information [27(. To 
derive equation of motion for p(t) we can use projection 
operator approach as following. We define a projection 
operator V by its action on an arbitrary operator A as 



PA = tr B {A}w, 



(13) 



where w is a pure bath operator satisfying the condition 
trs{w} = 1. For convenience we also define a comple- 
mentary projection operator Q = 1 — V . For the pro- 
jected density matrix VW(t) one can write a formally 
exact equation of motion as 

^VW{t) = -iVC x (t)PW(t) 
at 

-P£j(i)ex P j-i^ dTQd{T)Q^ QW(t ) 

- j dr'VLi(t)exp^-i £ dr"Q£/(r")e| 

x Q£j{t-T')rW{t-T'), (14) 

where the subscript j denotes interaction picture with re- 
spect to the bath Hamiltonian Hr. Eq. (fl"4"|) is known as 
Nakajima-Zwanzig identity [2^. l36l|. Although one can- 
not practically solve Eq. (fl4"|) , because it is as difficult as 
the original complete problem, one can use it as a starting 
point for approximations that will follow. 

B. Convolutionless equation of motion 



We start with the density operator W, which describes 
the state of the system composed of the aggregate and 
its surroundings. This operator fulfills Liouville-von Neu- 
mann equation 

^W{t) = - l -[H, W(t)}_ = -iCW(t). (10) 

Here, [A, £?]_ = AB — BA, and we introduced Liouville 
superoperator C by its action on an arbitrary operator 
A as LA = j-[H, Liouville superoperators £5, Lb 
and Ls-b that correspond to Hamiltonian operators of 



Before we introduce approximations into Eq. (|14[) . we 
can make one more formal step, which enables us to over- 
come its time non-local character. Let us assume an evo- 
lution superoperator W(t), which satisfies 

PW(t) =U'(t-t )PW(t ). (15) 

This superoperator is of course unknown, but it allows 
us to turn a time non-local differential equation (|14[) for- 
mally into a time-local one. We can write 

o 

—VW(t) = -iVd{t)VW{t) 



5 



-P£ z (f)exp j-ijf' drQ£ z (r)Q| CW(t ) 



t-to 



J dT'V£i(t)expl-i j dT"Qd(T")Q 



x Q£i(t - t')U '(-/) PW(t) 



(16) 



We do not develop the convolutionless (or time local) 
theory any further in this paper, because we will be in- 
terested only in terms up to the second order in £/. In 
such a case U'(—t) can be approximated as U'(—t) k, 
Us{— t)Ub{— t), where Us{—r) and Ub{—t) are evolu- 
tion superoperators with respect to Cq and £3. Inter- 
ested reader can refer e.g. to Refs. 25j [H, H(| for further 
details. 



C. System-Bath Coupling 

We will now assume the interaction Hamiltonian in a 
form of Eq. @ where index n now runs through all 
relevant single-exciton and multi-exciton states 



Hi = A ®nK n . 



(17) 



Correspondingly, K n = \u n ){u n \ for single excitonic 
states. Expanding Eqs. (fT4|) and (fTBj) into the second 
order in Hi we find that the third term on the right hand 
side (r. h. s.) can be conveniently expressed via so-called 
bath (or energy gap) correlation functions defined as 

C ron (r) =tr Q {U B (-T)A$ m U B (T)A$ n w eq }. (18) 

Here, we have chosen a specific form of the bath den- 
sity matrix w = w eq , where w eq is the canonical density 
matrix of the bath DOF. Defining also an operator 



A m (r) =J2C mn (r)U s (T)K n Ul(T) 

n 

and a superoperator MS 2 >(t) such that 

M^(r)A = Y}K m ,k m {T)Us{T)AU s {-T) 

m 

-Us(t)AU s (-t)AI(t)}-, 
the Eqs. (fT4")) and p^|) can be rewritten as 

d 



(19) 



(20) 



at 



P (t) = -iC sP {t) 



t-t 



]T / dr[K m ,A m (r)p(t)-p(t)AUr)U (21) 



and 



di 



p(t) = -iC s p(t) 



j dT{K m ,A m (T)U S (T)p(t~T)U S (-T) 



U s (T)p(t-T)Us(-T)Al(r)]-. 



(22) 



Provided we can supply a model for the correlation func- 
tion C mTl (r) we are in position to write down the equa- 
tions of motion for RDM in terms of known quantities. 
The last step necessary to implement these equations is 
to represent them in the basis of the eigenstates of the 
aggregate Hamiltonian. We define 



Pab(t) = (u a \p(t)\u b ), a,b=l,...,N, 



Pab(t) = (U a \p(t)\U b ), 



(23) 



a,b = N + 1,..., N + N(N - l)/2, (24) 

and in a similar manner for matrix elements of other 
operators and superoperators. This leads to 



0_ 

di: 



>ab(t) = -iuJ a bPab(t) - '^Habcdit) Pcd(t) , (25) 



with TZ a bcd(t) being the matrix elements of the superop- 
erator defined by the r. h. s. of Eq. (|2"Tj) . and 



di 



Pab(t) = -iu ab p ab {t) 



t-t 



cd 



M abcd {T)p cd {t - r), 



(26) 



with Ai a bcd(T) the matrix elements of the superoperator 
defined by the r. h. s. of Eq. (1221) . All the quanti- 
ties needed to calculate the matrix elements lZ a bcd(t) and 
M-abcd{t) are known provided the energy gap correlation 
function is know. 



D. Energy gap correlation function 

As a suitable model of the energy gap correlation func- 
tion we choose so-called multi-mode Brownian oscillator 
(BO) [3]. In general, Brownian oscillator model can inter- 
polate between underdamped intra-molecular DOF and 
(usually) overdamped bath DOF representing the imme- 
diate surroundings of the molecule. In this paper, we 
assume the correlation function of the energy gap of each 



6 



molecule in the aggregate to be the same, and indepen- 
dent of neighboring molecules, i.e. 



Cab(t) = C(t)S ab . 



(27) 



The correlation function C(t) is taken in a form of the 
overdamped BO model 

C(i) = -ih\ke~ m sgn t 



4AA ^ v n e 



-v n \t\ 



n—1 



A 2 



with 



= 2wn R 

Vn ~ hf3' p - k B T 



1 A 1 

, A = - 



(28) 



(29) 



Here, A is the reorganization energy, v n are so-called Mat- 
subara frequencies, fee is the Boltzmann constant, T is 
the thermodynamic temperature and r c is the so-called 
bath correlation time. The BO form of the correlation 
function satisfies all general constraints put of a correla- 
tion function by thermodynamics [§}. Apart of the tem- 
perature which we assume to be T — 300 K in all calcu- 
lations in this paper, the BO model is determined by two 
parameters only; by the reorganization energy A which is 
experimentally related to the Stokes shift S — 2A and by 
the bath correlation time r c . BO is a widely used, well 
physically motivated, but not the only possible model 
for the bath correlation function. Implications of other 
forms of the correlation function for the RDM dynamics 
will be studied elsewhere. 



E. Secular and constant relaxation rate 
approximations in the energy eigenstate basis 



Eqs. (|23|) and (|26p are systems of coupled (integro-) 
differential equations for the elements of the RDM. From 
the first terms on the r. h. s. we deduce that the ele- 
ment p a b{t) oscillates with a frequency close to w a f,. It 
is often justified to assume that two terms oscillating on 
different frequencies are independent of each other. For 
their envelopes p a b(t) = e lulabt p a b{t) we have 

^- t Pab(t) = -Y, Kabcdity^-^pcdit), (30) 



cd 



and integration over time has therefore a relatively 
smaller contribution when io a b — lu cc i ^ 0. Neglecting 
these contributions, usually termed secular approxima- 
tion [27J, leads to setting 



T^-abcd{t) — 0, 



(31) 



for all term except when a = c and b — d , or a = b 
and c = d. The interpretation of the remaining non-zero 



terms is simple. The terms lZ aa b b (t) for a ^ b repre- 
sent rates of transition from level denoted by index b to 
a level denoted by a. The term lZ aaaa (t) corresponds to 
the total transition rate from the level a to all other lev- 
els. The terms TZ a bab(t) (a ^ b) are rates of the damping 
of a coherence element p a \, (t) . In secular approximation, 
the dynamics of populations of electronic levels is thus 
decoupled from the dephasing of coherences. The above 
arguments for the secular approximation apply also to 
the integro-differential equation (f2l)|) . and we can thus 
define four different second order equations of motion for 
the RDM, with different levels of approximation. From 
the perspective of our derivation, the most general sec- 
ond order equation is Eq. (|2l))) . which we have denoted 
full TNL. The convolutionless Eq. flU]) denoted full TL 
can be regarded its approximation, but it can also be 
alternatively viewed as derived by different cummulant 
approximation, see Refs. [25|,[2(|. The set of four meth- 
ods investigated here is completed by applying secular 
approximation to the full TNL and full TL equations. 

All four sets of equations of motion we consider here 
are extensions to the two well-known constant relaxation 
rate theories. To arrive at the well-known Redfield equa- 
tions [22}, one can assume certain coarse graining of the 
RDM dynamics so that all significant changes to the p(t) 
occur on a time scale much longer than the correlation 
time r c . Then time to in Eq. (|2ip can be put to —00 and 
the integration limits are then from zero to infinity. The 
relaxation tensor 1Z thus becomes time independent. If 
we, on the other hand, consider the decay of C(t) to be 
much faster than even the transition frequencies between 
electronic levels, we can assume C(t) w CoS(t) and Eq. 
(|2"Tj) has the well-known Lindblad form [23, |37|. Only 
for the Lindblad form and for the Redfield equations in 
secular approximation, it can be shown that the diago- 
nal elements of p(t) are always positive. For all other 
equations we have derived here, this assertion cannot be 
proven in general. This is a consequence of the fact that 
they are derived in a low order of perturbation theory. 



IV. NON-LINEAR SPECTROSCOPIC SIGNALS 

Non-linear spectroscopic signals are very well described 
by time-dependent perturbation theory [3]. The equa- 
tions of motion, Eqs. (|2"Tj) to (|2"2")l . can be extended by 
semiclassical light-matter interaction term. This yields 



d_ 

dt 



p(t) = -iC sP (t) - V[p{t)]{t) + iVp(t)E(t), (32) 



where E(t) = n ■ E(t), is the projection of the external 
electric field vector E(t) on the normal vector n in di- 
rection of the molecular transition dipole moment. The 
symbol D[p(t)](t) represents the relaxation term chosen 
from the full TNL, full TL, secular TNL and secular TL 
equations of motion. The superoperator V is a commu- 
tator with the dipole moment operator fj, — np, so that 
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for an arbitrary operator A we have 

VA=~\p,A]-. (33) 
a 

A. Third-order non-linear response theory 

Non-linear optical signals are related to the RDM via 
polarization 

P(t) = tr{w{t)}. (34) 

In particular, for the third order non-linear signal E^ (t) 
one can write 

Ei 3 \t) » ^^P (3) (^) = iu *r{/*p (3) (t)}, (35) 

where the upper index ( 3 ) denotes that the quantity is of 
the third order of the perturbation theory with respect to 
the external electric field E(t). By defining the evolution 
superoperator U{t) which fulfills Eq. with E{t) = 
we can write the third order perturbation term as 

oo oo oo 

P m (t) = -l j J J dT 3 dT 2 dTlU(T 3 )VU(T 2 )VU(n)V P 


X E(t - T 3 )E{t - T 3 - T 2 )E{t -T3-T2- n). (36) 

In experiment, the laser field is often prepared in a form 
of three incident pulses 

E(t) =A x {t- tl ) e - ini{t - tl)+lkir 
+A 2 {t - t 2 )e- in * {t - t2)+lk2r 



t) we can write for the time and the delay dependent 
signal field 

E s (t,T,r) « R 2g (t,T,r) 
+ R 3g (t, T, t) - R* lf (t, T, t), t ^ 0, (38) 
E s (t,T,T)*R lg (t,T,\r\) 

+ R ig (t,T,\r\) - R% f (t,T,\T\), r < 0. (39) 

The absolute value in Eq. ([3"9"[) originates from the fact 
that response functions R are defined for positive time 
arguments only, and negative r is achieved by switching 
the order of the k\ and k 2 pulses. The individual re- 
sponse functions R are listed in the Appendix |XJ Most 
importantly, they consist of series of propagation of the 
density matrix blocks by evolution operators obtained 
from the solution of equations of motion. We have e.g. 

R 2g (t,T, T ) =tr{fx ge U egeg (t)Vi^ 

x U eeee (T)V e ^U gege (r)V g ^p }, (40) 

where the evolution superoperators act on an arbi- 
trary operator A as a dipole operator fi a f, from the right, 
i.e. V^f* A = Apai,. The superoperator V$ is defined 
analogically with the action of p a b from the left. The in- 
dices e and g denote electronic bands as denoted in Fig. 
[TJ Thus, the above operators and the action of superop- 
erators on an arbitrary operator A are expressed in the 
basis of Hamiltonian eigenstates as 

po = \g)(g\, (41) 



+ A 3 (t ~ t a ) e -<na(t-t*)+a*r + c c ^ ( 37 ) 

with different k- vectors fci, k 2 and £3. In the rest of the 
paper we assume Q\ — Q 2 — Q 3 = O, A\{t) = A 2 (t) = 
A 3 (t) = A(t). The expression obtained by inserting Eq. 
([3^)1 into Eq. ([33)1 can be significantly simplified in cases 
where the system consists of a ground-state and a band of 
excited states, with the transition frequency close to reso- 
nance with the laser pulse frequency fi, and by assuming 
the laser pulses are ultra short, i.e. A(t) w Eo5(t). For 
an experiment which detects non-linear signal emitted in 
the direction — fei + k 2 + k 3 , the third order signal has 
a frequency w f2 and it is obtained from just a handful 
of response functions that represent certain contributions 
to the triple commutator in Eq. )[3l)|). The details of the 
derivation can be obtained e.g. in Ref. fTTj | . 

If the delays between the pulses are selected such that 
r denotes the delay between the first (ki ) and the second 
(£2) pulses, and T denotes the delay between the second 
and third (k 3 ) pulse (e.g. £3 = 0, t 2 = —T and t\ = — T — 



^e g =J2fn g 9) K)(9\, (42) 

71 

U gege {t)[A] = ]T U^(t){u m \A\g)\g){u n \, (43) 



u eeee {t)\A]= £ ui°zl(t) 

nn'mm' 

x {u m >\A\u m )\u n ){u n >\- (44) 

Eqs. (|4"Tj) to (|4"4"]l together with the App endix lAl enable us 
to calculate expected non-linear signal from the knowl- 
edge of the matrix elements of the evolution superop- 
erator. This type of knowledge can be obtained from 
solutions of the four different equations of motion that 
we presented in Section [Ml 
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Table I: Parameters of the model trimer. The parameter e n 
represents the transition energy of n— th monomer, transition 
dipole moments d n are taken relative to some value do- Pa- 
rameters h n and a n axe explained on Fig. [2] 

B. Two-dimensional Coherent Spectroscopy 

Two-dimensional coherent spectrum, S(wt,T, ui T ), is 
obtain from the non-linear signal by Fourier transforming 
the time and pulse delay dependent signal electric field 
Es(t, T, t) along the t and r variables [1, [ll[ as 

oo oo 

E(u u T,u T )= J dt J dTE s {t,T,T)e iUtt - iu ^. (45) 

—00 — 00 

The Fourier transform in r yields an uj t dependence that 
is formally similar to linear absorption spectrum, while 
the transform in t yields generalized absorption and stim- 
ulated emission from a non-equilibrium state created by 
the first two laser pulses. 2D spectrum thus represents a 
2D absorption/emission and absorption/absorption cor- 
relation plot. During the pulse delay time T the system 
evolves both in the electronically excited state and in the 
ground state, but no optical signal is generated. Relax- 
ation of populations in the electronically excited band 
leads to evolution of non-diagonal 2D spectral features, 
so-called cross-peaks. Cross-peaks appearing at T = 
are a signature of excitonic origin of the observed ex- 
cited states. The 2D cross-peaks oscillate in T as long 
as the corresponding electronic coherence elements of the 
reduced density matrix are oscillating. The life time of 
the electronic coherences can thus be estimated directly 
from the T dependent sequence of 2D spectra [IH, [TJ • 

V. NUMERICAL RESULTS AND DISCUSSION 

In this section we study dynamics of model aggregate 
viewed via population and coherence dynamics and via 
2D coherent spectrum. We define a simple model ag- 
gregate for which we calculate excited state dynamics in- 
cluding evolution of coherences between electronic states, 
linear absorption and 2D spectra at chosen population 
times. Calculations of linear absorption, which require 
only knowledge of the time evolution of optical coher- 
ences, are performed using the secular time local equa- 
tion, since it is known to yield exact result at least for 
some models j35|. Population dynamics is calculated us- 
ing all four methods we discussed in Section llHl and the 
results are compared. 

The simplest model of an aggregate that can exhibit 
all effects observed in Ref. [1J| is a trimer. The ge- 



ometry of the models, together with the meaning of the 
parameters is presented in Fig. [2l In Tab. Uwe summa- 
rize the main parameters of the model. Parameters J, h 
and do from Tab. |T] are not independent. For given h n 
and d n we could in principle calculate the value of reso- 
nance coupling J. Because we are not interested in the 
absolute amplitude of the absorption or 2D spectra we 
assume do to be fixed by the values of h n and d n to yield 
the expected value of J. All three resonance couplings 
J between the molecules are set to J = 200 cm -1 for 
the calculations presented here. The values of the tran- 
sition dipole moments determine the initial condition for 
the population dynamics. We assume that the excita- 
tion light intensity and the value of the transition dipole 
moment are such that the system is only weakly excited. 
The total population of the excited state band is normal- 
ized to 0.01. The relative values of the transition dipole 
moments are chosen so that the linear absorption spec- 
trum (see Fig. \3§ shows peaks of roughly the same height. 
Two peaks originating from the energetically lowest and 
the energetically highest states dominate the spectrum, 
the third level contributes as a shoulder to lowest energy 
peak. 

Two parameters that influence the coupling for the 
model system to the bath are reorganization energy A 
and correlation time r c . We vary these parameters in the 
range that can conceivably represent chlorophylls in pho- 
tosynthetic complexes (see e. g. Refs. (38l l39l|). 



A. Population relaxation and evolution of 
coherences 

First, we compare relaxation dynamics of populations 
of excited state of our aggregate after excitation by an 
ultrashort laser pulse. TL equations of motion where 
solved by standard numerical methods for ordinary dif- 
ferential equations provided by the Mathematica® soft- 
ware. For the TNL equations we used fast Fourier trans- 
form method. Figure @] presents the first 1 ps of the 
population dynamics after a 5— pulse excitation of the 
trimer from Tab. [I] at the temperature T = 300 K. Re- 
organization energy A = 120 cm -1 and correlation time 
t c = 50 fs are the same at all three monomers. The 
dynamics with the same parameters for a selected coher- 
ences element pia(t) is presented in Fig. [5] The overall 
conclusion is that all four methods yield a similar general 
behavior for the populations, with some difference at the 
short time evolution and also slightly different long time 
equilibrium. Examination of the Figure pleads us to the 
conclusion that the methods yield two different results - a 
short coherence life time for the time local methods, and 
a relatively longer life time in case of the time non-local 
methods. The behavior of the coherence pi3(t) repre- 
sents a general tendency that we have observed for all 
electronic coherences over a wide range of parameters. 

Let us now concentrate on short time behavior of the 
populations and coherences in more retail. In the short 
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Figure 2: Geometry and parameters of a trimer aggregate. 
One monomer is chosen to be positioned at the origin of the 
coordinate system, with the transition dipole moment point- 
ing along the x axis. The positions the transition dipole mo- 
ments of the other two molecules in space are characterized 
by their distance hi and /13 from the origin of coordinates 
and by the angles 0:2 and 03 . Orientations and lengths of the 
dipoles are given in Tab. [I] In our example we assume that 
the aggregate is planar. 
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Figure 3: Linear absorption spectrum of the model trimer 
for various parameters of the system bath interaction: (a) 
A = 120 cm -1 , t c = 50 fs, (b) A = 30 cm -1 , r c = 100 fs, 
calculated by the secular TL theory (full lines) and the secular 
TNL theory (dashed lines). 



time evolution of the coherences the four methods group 
into two distinct groups with short (TL methods) and 
long (TNL methods) coherence life time. Whether the 
underlying equation is secular or not seems to have only a 
little influence on the coherence dynamics. Fig. [5] shows 
the short time (0 — 400 fs) comparison of the popula- 



tion calculated by the four equations of motion. We can 
clearly see that the results can be naturally grouped ac- 
cording to the presence of fast oscillatory modulation of 
the population relaxation dynamics. In the one group 
we have the full TL and full TNL methods, where such 
oscillations clearly occur, the second group comprises the 
two secular methods with no oscillations present. Thus, 
it can be concluded that the non-secular terms in the 
equations of motion are the cause of these oscillations. 
This is also supported by comparison of the population 
dynamics of the full TL and full TNL equations from Fig. 
IU(e.g. the population of the state 1). The oscillation on 
the full TNL curve last longer than those of the full TL 
one, which reflects the longer coherence life time we have 
found for the TNL equations. 

Let us now discuss the long time limit of the time evo- 
lution. As expected, the two secular theories yield the 
same equilibrium at long population times. This equilib- 
rium corresponds to the canonical distribution of popu- 
lation among the excitonic levels at T = 300 K. In both 
secular TNL and secular TL cases, coherences have re- 
laxed to zero at long times as the inset of the Fig. [3] 
demonstrates. The non-secular TNL and TL equations 
yield non-zero, stationary coherences at long times, and 
correspondingly, the long time equilibrium populations 
do not correspond to the canonical thermal equilibrium. 
Although both non-secular theories converge to results 
different from the canonical equilibrium, the full TNL 
equation yields populations that are physical at all times 
for the studied system parameters, i.e. they are always 
positive. The full TL equation on the other hand fails to 
keep probabilities positive at long times, and the occu- 
pation probability of the highest electronic level becomes 
negative after 200 fs for the parameters used on Fig. |4j 

In light of recent experiments [lij , the conclusion that 
time non-local theories lead to a longer coherence life 
time than the time-local ones (i.e. also longer than the 
standard constant rate theories) is probably the most in- 
teresting. We have performed calculations of the RDM 
dynamics while varying the reorganization energy and 
the correlation time. The absolute values of the coher- 
ence /0i3 (t) elements were fitted by a single exponential 
to estimate coherence life-time. The results are summa- 
rized in Fig. [7J The Fig. [JJ\. shows the results for secular 
TL and secular TNL equations. Clearly, with growing 
correlation time r c , the full TNL equations lead to a in- 
creasing coherence life time. The full TL equation shows 
only a very weak dependence of the coherence life time 
on correlation time. Another interesting observation is 
that for correlation time longer then 50 fs, the depen- 
dence of the coherence life time on the reorganization 
energy A is different for full TNL and TL methods. Time 
local theory, in accordance with the standard rate the- 
ories, predicts decrease of the coherence life time with 
A. The full TNL theory predicts (within the parameter 
range studied here) an opposite tendency. The Fig. [7)3 
shows similar conclusion for the non-secular versions of 
the theories, with the same difference between TL and 
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Figure 4: First 1000 fs of the excited state population dy- 
namics of a trimer with parameters A = 120 cm -1 , r c = 50 
fs, calculated by all four methods. For these particular pa- 
rameters, the full TL equation breaks positivity of the RDM 
diagonal elements after 200 fs. Its prediction for the popula- 
tions of the lowest and highest levels is significantly different 
from the other three methods. 
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Figure 5: First 500 fs of the dynamics of the RDM coherence 
element pi3(t), with parameters from Fig. [4] calculated by 
all four methods. Detail of the long time part of the time 
evolution is presented in the inset. 



TNL theory. The dependence of the coherence life time 
on A in case of TNL equations is not monotonous. 



B. Two-dimensional spectrum 

As discussed in the Introduction, the secular TL equa- 
tion of motion yields an exact result for the dephasing 
of an isolated optical coherence [35]. One can show, by 
comparison of the absorption spectra calculated by secu- 
lar TL and TNL methods (see Fig. ©, that the TNL 
theory leads to certain artifacts (second peak) and is 
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Figure 6: First 400 fs of the population dynamics of the trimer 
with parameters A = 30 fs and r c = 100 fs. Results of full TL 
and TNL theories are presented in upper subfigure (A), the 
secular results are found the the lower subfigure (B). 



therefore not suitable for the description of the optical 
coherence evolution. Consequently, on can only hope to 
obtain valid results for the evolution superoperators at 
the first and the third time interval of the third order re- 
sponse functions by the TL theories. In Ref. [40j it was 
shown that non-secular terms in the TL equations for op- 
tical coherences lead to temperature dependence of the 
positions of excitonic bands in absorption spectra. This 
dependence was shown to be strong when the electronic 
states involved are characterized by significantly different 
reorganization energy [34|,[4(|. Indeed it can be shown for 
homodimer that the non-secular terms are exactly zero 
in second order TL theory if the monomers exhibit the 
same reorganization energy [34J. We can therefore ex- 
pect the non-secular effects in the optical coherences to 
be weak in our case, and we choose secular TL to cal- 
culate the evolution superoperators in the first and third 
time interval of the response function, Eq. (|40[) . 

Concerning the population interval, the situation is 
somewhat different. As we have shown above, the non- 
secular TL theory leads to dynamics that breaks the pos- 
itivity condition for the population probabilities at long 
times. At the same time, short time dynamics is very sim- 
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Figure 7: The life time of coherence pi3{i) as obtained from 
fitting the coherence dynamics calculated by all four methods 
for various parameters A and r c . The upper subfigure (A) 
shows the life times obtained by the secular methods, while 
the lower subfigure (B) presents the same for non-secular 
methods. 



ilar to the full TNL. Both theories predict population os- 
cillations during the life time of the electronic coherences. 
The full TNL equation, however, preserves positivity, at 
least for the parameters studied here, and can be there- 
fore used to calculate meaningful 2D spectra. For the 
same reason, both secular theories can also be success- 
fully used to calculate 2D spectrum. As the oscillation 
of the populations predicted by non-secular theories are 
too small to be reliably observed in 2D spectrum (only a 
small change of the crosspeak amplitude due to the pop- 
ulation transfer is observed after 140 fs of relaxation in 
2D spectrum of Fig. [5]) we expect only a small difference 
of the 2D spectrum to appear between the secular and 
full TNL theories. For the calculation of the representa- 
tive 2D spectrum we therefore choose the secular TL and 
the full TNL theories. These two differ from each other 
mainly in the prediction of the life time of the electronic 
coherences. The observable difference in the calculated 
2D spectra should therefore predominantly result from 
the different life time of the electronic coherence. 
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Figure 8: Two-dimensional coherent spectra of the trimer 
model at population times T = 0, 20 and 140 fs calculated 
by the secular TL method (left column) and the full TNL 
method (right column). The system-bath interaction param- 
eters are A = 30 cm -1 and r c = 100 fs. The coherence element 
piz{t), which is mainly responsible for the oscillatory behav- 
ior of the crosspeaks, is presented in the upper right corner 
of the figure. The 2D spectrum at T = fs is the same 
for both methods and is therefore presented only once. The 
population times are selected so that they represent different 
phases of the pi3(i) element (denoted by arrows on the co- 
herence element figure). Arrows in the 2D spectra denote the 
orientation of the peaks. All spectra are normalized to 1 with 
contour step of 10 %. Positive features are in full red line, 
negative features are represented by dashed blue line, and the 
zero contour is depicted by the full black line. 



Fig. [8] presents 2D spectra for A = 30 cm -1 and 
t c = 100 fs. These parameters lead to a rather slow 
relaxation and consequently to narrow spectral peaks in 
both absorption (see Fig. [3]) and 2D spectra. This allows 
us to clearly see characteristic T— dependent oscillations 
of the peaks in 2D spectrum. At T — fs, both methods 
provide the same 2D spectrum, with four peaks. Two 
diagonal peaks arise when all three perturbations of the 
system by electric field occur on the same level, while two 
crosspeaks appear from interactions occurring on differ- 
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ent levels. Negative peaks correspond to excited state 
absorption (see Fig. [!}. F° r two molecules that are not 
excitonically coupled, all contributions to the crosspeaks 
cancel out exactly, while if two molecules are excitonically 
coupled non-zero crosspeaks appear. The shapes of the 
peaks are influenced by the phase evolution of the coher- 
ence elements of RDM during the population time T. On 
the upper left figure of Fig. [5] we have marked the elon- 
gation of the diagonal and off-diagonal peaks by arrows. 
The elongation can be best judged by looking at the zero 
contour (in black). This particular elongation is char- 
acteristic for the phase of the pi3(t) element (see upper 
right figure of Fig. |S) at T = 0. At T = 20 fs the phases 
of the pisit) calculated by both methods are opposite to 
the phase at T — 0. The 2D spectra calculated by the two 
different methods at T — 20 fs differ only in the precise 
positions of the contours. This phase of the coherence 
element is characterized in 2D spectrum by a different 
orientation of the peaks. Interestingly, at T = 140 fs the 
two methods predict pisit) that have mutually opposite 
phases and as a consequence the 2D spectra at T = 140 fs 
calculated by different methods differ in the orientation 
of their crosspeaks. Since the secular TL theory pre- 
dicts a simple dephasing of the coherence and a regular 
oscillation with a single frequency proportional to the en- 
ergy difference between corresponding energy levels, it is 
in principle possible to distinguish, even experimentally, 
deviations from this prediction. Our conclusion is that 
such a deviation should be a consequence of the memory 
effects in the reduced system time evolution. 



C. Validity of secular and Markov approximations 

Several conclusions about the applicability of the sec- 
ular and Markov approximations can be drawn from the 
above results. As pointed out in Ref. [3^, Markov ap- 
proximation, which in the second order in system-bath 
coupling converts the TNL equations to the TL ones, 
leads accidentally to an exact result for an optical co- 
herence element interacting with the harmonic bath. It 
has been also pointed out previously [33l |4l| that in the 
same case, the TNL equations lead to artifacts. When 
studying relaxation dynamics of the populations and elec- 
tronic coherences in excitonic systems, full TL theory 
leads to a breakdown of the positivity of the RDM, while 
none of the secular theories suffer from this problem. In 
principle, the full TNL theory suffers from this problem, 
too. However, it has been found less susceptible to it 
here. The secular theories lead to canonical density ma- 
trix at long times, while the full TNL results in a sta- 
tionary state characterized by non-zero (but constant) 
coherences. Such result corresponds to an additional 
renormalization of the electronic states by the interac- 
tion with bath , an d has to be expected even at a weak 
coupling limit |42j. It is important to note in this con- 
text that the canonical equilibrium is to be expected for 
the system consisting of the molecule and the bath as a 



whole, not for its parts [42| . 

For the population dynamics we are therefore forced 
to conclude that the full TNL theory represents the best 
candidate for a correct description of relaxation phenom- 
ena in the second order of the system bath interaction. It 
predicts similar population transfer times as other meth- 
ods, it is much less sensitive to the breakdown of the 
positivity than its TL counterpart, and it leads to a 
bath renormalization of the canonical equilibrium. Most 
interestingly however, it predicts longer coherence life 
time than the TL theory. It was recently established 
by Ishizaki and Fleming [33| that this is to be expected 
from a higher order theory. 

In the light of the above conclusions about the dy- 
namics of optical coherences and the populations and co- 
herences of the one exciton band, we suggest a hybrid 
approach to calculating 2D spectra, which consists of the 
application of the TL method on optical coherences (first 
and third time interval) and the full TNL method on the 
calculation of the RDM dynamics in the one exciton band 
during the population time T. 



VI. CONCLUSIONS 

In this paper we have compared four different theories 
of excitation energy transfer and relaxation in molecu- 
lar aggregate systems, with a special attention paid to 
lifetime of electronic coherences. Second order time non- 
local and time local theories with and without secular 
approximation were studied. For our specific model of 
an aggregate we have concluded that time non-local the- 
ories can account for experimentally observed electronic 
coherence life time that is significantly longer than the 
one predicted by the standard time-local secular relax- 
ation rate theory. Markov approximation leading to time 
local equations of motion was found to be responsible for 
the reduction of the coherence lifetime, while the influ- 
ence of the secular approximation on the life time was 
found rather weak. The time local theory without secu- 
lar approximation is found to break positivity of the oc- 
cupation probabilities in the range of parameters studied 
here. We conclude that time-local second order theory 
is not suitable for simulating the coherence transfer ef- 
fects. Simulations of two-dimensional spectra show that 
the time non-local effects can be experimentally identi- 
fied based on the analysis of the oscillations of the cross 
peaks. 
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Appendix A: Third Order Response Functions 

In this appendix we list the third order response func- 
tion used in calculating the impulsive 2D spectra. The 
first index of the response function follows the standard 
notation of Ref . [3] . The second index is g for pathways 
not involving the two-exciton band, while all pathways 
denoted by f include a two-exciton contribution (see e.g. 
Ref. [HI). 

R lg (t,T,r) =tr{p ge U eg e g (t)V^ 



x U gggg {T)VifU gege {T)VWpo}, (A3) 



R 4g (t,T, T ) = tr{fi ge Ue geg (t)V^ 



x U gagg {T)V<pU egBg {T)VWp }, (A4) 



R lf (t,T,T) = tr{ f Xf e U efe f(t)vff 



x U eeee (T)V$U egeg (T)V&po}, (Al) 



R 2g (t,T,r) =tr{fi ge Ue geg (t)vW 



x U eeee (T)V^U egeg (r)V^ Po } (A5) 
R 2f (t,T,T) =tr{p fe U ef ef(t)V { e f 



x U eeee (T)V^U gege (r)V^ Po }, (A2) 



R 3g (t,T,r) = tr{p ge U egeg (t)V^ 



x U eeee (T)V^U gege (r)V^p }. (A6) 

Operators and superoperators used in this appendix are 
defined in Section [IV Al 
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